#=
对B^{-K}中的项目求导
原则上src下第一级别的代码之间互相不依赖，adso2.jl可以理解成asdo.jl的后半部分
=#

"""
B矩阵对mu的导数
"""
function pbpmu_calc_tapemap(sp::Splitting, mui::Vector{Float64})
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 2*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            #求出没有mu的时候的ehk
            #后面作用mu上去，需要在没有mu的ehk上作用，才体现的是加了mu以后B矩阵的变化，否则就是已经有mu的B又加了mu
            ehkslice = exp(sp.dt*mui[xi])*ehk2[[xi, xi+Nx], :]
            #如果认为mu出现在右侧的话，那么要修改的就是列，比较麻烦
            funcforward1 = (mu) -> bmat_IsingADX(exp(-sp.dt*mu)*ehkslice, Val(-1), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            funcforward2 = (mu) -> bmat_IsingADX(exp(-sp.dt*mu)*ehkslice, Val(1), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            jacxi1 = derivative(funcforward1, mui[xi])
            jacxi2 = derivative(funcforward2, mui[xi])
            tapemap[xi, 2*si-1] = AutoJacTape([mui[xi]], jacxi1; val=-1, mu=mui[xi])
            tapemap[xi, 2*si] = AutoJacTape([mui[xi]], jacxi2; val=1, mu=mui[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵对格点i上面的mu的导数
"""
function pbpmu_calc_xi(sp::Splitting, hsxi::Vector{Int}, xidx::Int, mui::Vector{Float64}; tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpmu2 = zeros(Nt, 2, 2*Nx)
    pbipmu2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = hsxi[bi] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 4, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpmu2[bi, :, :] = jac1[1:2, :]
            pbipmu2[bi, :, :] = jac1[3:4, :]
        end
    else
        @warn "计算jacobian"
        #
        for bi in Base.OneTo(Nt)
            ehk2, Ui = unpack_splitting(sp, bi)
            ehk2 = Diagonal(exp(sp.dt*mui))*ehk2
            #因为bmat_IsingADX中不接受vector的输入，jacibian不会重新定义exp（只有exp.）
            #这里要用derivative
            funcforward2 = (mu) -> bmat_IsingADX(exp(-sp.dt*mu)*ehk2[[xidx,xidx+Nx], :],
            Val(hsxi[bi]), Ui[xidx], 0.0, 0.0, 1.0, sp.dt)
            #
            jacxit = derivative(funcforward2, mui[xidx])
            jac1 = reshape(jacxit, 4, 2*Nx)
            pbrpmu2[bi, :, :] = jac1[1:2, :]
            pbipmu2[bi, :, :] = jac1[3:4, :]
        end
    end
    return pbrpmu2, pbipmu2
end


"""
计算∂L/∂mu
"""
function meas_gradmu(ss::ScrollSVD{T}, allbmats, sp, allcfgs, mus; tapemap=missing) where T
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    #这里必须指定tapemap，且必须是用baress算出来的
    if ismissing(tapemap)
        throw(DomainError("tapemap"))
    end
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    muibar = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        rpmu, ipmu = pbpmu_calc_xi(
            sp, allcfgs[:, xi], xi, mus; tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            muibar[xi] += sum(rpmu[bi, 1, :] .* real(bbar[bi, xi, :]))
            muibar[xi] += sum(rpmu[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            muibar[xi] -= sum(ipmu[bi, 1, :] .* imag(bbar[bi, xi, :]))
            muibar[xi] -= sum(ipmu[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
    end
    return muibar
end


#=
处理tp的方法1:在右侧乘exp(dt*dp)
缺点是每一个tp都会修改所有的ehk[i,:]，因为乘在右侧修改的是列
ehkslice = ehkbare*(I + offdiag)[i,:] = ehkbare[i,:] + ehkbare*offdiag[i,:]
如果offdiag链接的是l,m两个格子
ehkbare[i,l] = ehkbare[i,m]*offdiag[2,1] + ehkbare[i,l]*offdiag[1,1]
ehkbare[i,m] = ehkbare[i,l]*offdiag[1,2] + ehkbare[i,m]*offdiag[2,2]
需要做的就是在ehkslice上面加入
exp([0 dt*tp; dt*tp 0]) = [cosh(dt*tp) sinh(dt*tp); sinh(dt*tp) cosh(dt*tp)]
function run()
ehkslice = rand(2,6)
func = function (x)
arr = hcat([x+cos(x)],zeros(1,4),[x+sin(x)]) 
arr = vcat(arr, arr)
return ehkslice + arr
end 
derivative(func, 1.0)
end
=#
#=
处理tp的方法2:在左侧乘exp(dt*dp)
缺点是每一个tp要修改两行的ehk
ehkslice = ((I + offdiag)*ehkbare)[i,:] = ehkbare[i,:] + (offdiag*ehkbare)[i,:]
如果offdiag链接的是i,j两个格子
ehkbare[i,:] = offdiag[1,2]*ehkbare[j,:] + offdiag[1,1]*ehkbare[i,:]
ehkbare[j,:] = offdiag[2,1]*ehkbare[i,:] + offdiag[2,2]*ehkbare[j,:]
=#
#=
这里采用第二种方案，只计算某两个格点的tp在所有ehk上产生的效果
=#

"""
B矩阵对hopping矩阵中非对角项的导数
"""
function pbptp_calc_tapemap(sp::Splitting, offidx::Vector{Int}, offval::Vector{Float64})
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 2*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            i1 = xi
            i2 = offidx[i1]
            ehkslice1 = ehk2[[xi, xi+Nx], :]
            ehkslice2 = ehk2[[i2, i2+Nx], :]
            #左侧乘一个exp（-dt*tp）逆矩阵
            #跟mu的情况类似，需要加上tp之后的变化，而ehk中已经有了，就需要暂时先去掉tp
            ehkslice1 = cosh(sp.dt*offval[xi])*ehkslice1 + sinh(sp.dt*offval[xi])*ehkslice2
            ehkslice2 = cosh(sp.dt*offval[xi])*ehkslice2 + sinh(sp.dt*offval[xi])*ehkslice1
            #
            funcforward1 = (mu) -> bmat_IsingADX(ehkslice1*cosh(-sp.dt*mu)+ehkslice2*sinh(-sp.dt*mu), Val(-1), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            funcforward2 = (mu) -> bmat_IsingADX(ehkslice2*cosh(-sp.dt*mu)+ehkslice1*sinh(-sp.dt*mu), Val(-1), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            funcforward3 = (mu) -> bmat_IsingADX(ehkslice1*cosh(-sp.dt*mu)+ehkslice2*sinh(-sp.dt*mu), Val(1), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            funcforward4 = (mu) -> bmat_IsingADX(ehkslice2*cosh(-sp.dt*mu)+ehkslice1*sinh(-sp.dt*mu), Val(1), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            jacxi1 = derivative(funcforward1, offval[xi])
            jacxi2 = derivative(funcforward2, offval[xi])
            jacxi3 = derivative(funcforward3, offval[xi])
            jacxi4 = derivative(funcforward4, offval[xi])
            tapemap[xi, 2*si-1] = AutoJacTape([offval[xi]], vcat(jacxi1, jacxi2); val=-1, offval=offval[xi])
            tapemap[xi, 2*si] = AutoJacTape([offval[xi]], vcat(jacxi3, jacxi4); val=1, offval=offval[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵对格点i上面的mu的导数
"""
function pbptp_calc_xi(sp::Splitting, hscfg::Matrix{Int}, xidx::Int, offidx::Vector{Int}, offval::Vector{Float64};
    tapemap=missing)
    #
    Nt, Nx = size(hscfg)
    pbrpmu2 = zeros(Nt, 4, 2*Nx)
    pbipmu2 = zeros(Nt, 4, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = hscfg[bi, xidx] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 8, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpmu2[bi, 1:2, :] = jac1[1:2, :]
            pbipmu2[bi, 1:2, :] = jac1[3:4, :]
            #另外一个位置
            tape = hscfg[bi, offidx[xidx]] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 8, 2*Nx)
            pbrpmu2[bi, 3:4, :] = jac1[5:6, :]
            pbipmu2[bi, 3:4, :] = jac1[7:8, :]
        end
    else
        @warn "计算jacobian"
        #
        for bi in Base.OneTo(Nt)
            ehk2, Ui = unpack_splitting(sp, bi)
            i1 = xidx
            i2 = offidx[xidx]
            ehkslice1 = ehk2[[xidx, xidx+Nx], :]
            ehkslice2 = ehk2[[i2, i2+Nx], :]
            ehkslice1 = cosh(sp.dt*mu)*ehkslice1 + sinh(sp.dt*mu)*ehkslice2
            ehkslice2 = cosh(sp.dt*mu)*ehkslice2 + sinh(sp.dt*mu)*ehkslice1
            #因为bmat_IsingADX中不接受vector的输入，jacibian不会重新定义exp（只有exp.）
            #这里要用derivative
            funcforward1 = (mu) -> bmat_IsingADX(ehkslice1*cosh(-sp.dt*mu)+ehkslice2*sinh(-sp.dt*mu), 
            Val(hsxi[bi]), Ui[xi], 0.0, 0.0, 1.0, sp.dt)
            funcforward2 = (mu) -> bmat_IsingADX(ehkslice2*cosh(-sp.dt*mu)+ehkslice1*sinh(-sp.dt*mu),
            Val(hsxi[bi]), Ui[xidx], 0.0, 0.0, 1.0, sp.dt)
            #
            jacxit1 = derivative(funcforward1, offval[xidx])
            jacxit2 = derivative(funcforward2, offval[xidx])
            jac1 = reshape(jacxit1, 4, 2*Nx)
            jac2 = reshape(jacxit2, 4, 2*Nx)
            pbrpmu2[bi, 1:2, :] = jac1[1:2, :]
            pbrpmu2[bi, 3:4, :] = jac2[1:2, :]
            pbipmu2[bi, 1:2, :] = jac1[3:4, :]
            pbipmu2[bi, 3:4, :] = jac2[3:4, :]
        end
    end
    return pbrpmu2, pbipmu2
end


"""
计算∂L/∂tp
"""
function meas_gradtp(ss::ScrollSVD{T}, allbmats, sp, allcfgs, offidx, offval; tapemap=missing) where T
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    #这里必须指定tapemap，且必须是用baress算出来的
    if ismissing(tapemap)
        throw(DomainError("tapemap"))
    end
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    muibar = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        rpmu, ipmu = pbptp_calc_xi(
            sp, allcfgs, xi, offidx, offval; tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            muibar[xi] += sum(rpmu[bi, 1, :] .* real(bbar[bi, xi, :]))
            muibar[xi] += sum(rpmu[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            muibar[xi] += sum(rpmu[bi, 3, :] .* real(bbar[bi, offidx[xi], :]))
            muibar[xi] += sum(rpmu[bi, 4, :] .* real(bbar[bi, offidx[xi]+Nx, :]))
            muibar[xi] -= sum(ipmu[bi, 1, :] .* imag(bbar[bi, xi, :]))
            muibar[xi] -= sum(ipmu[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            muibar[xi] -= sum(ipmu[bi, 3, :] .* imag(bbar[bi, offidx[xi], :]))
            muibar[xi] -= sum(ipmu[bi, 4, :] .* imag(bbar[bi, offidx[xi]+Nx, :]))
        end
    end
    return muibar
end


"""
计算包含额外权重的导数
``````
dlnw = dlnη + dlndet(I+B) = d∑_{τ,i}ln η(τ,i) + dlndet(I+B)

dlnw / dθ_{j} = d∑_{τ,i}ln η(τ,i)/dθ_{i} + dlndet(I+B)/dθ_{j}
= d∑_{τ}ln η(τ,j)/dθ_{j} + ∑_{τ}dlndet(I+B)/dB^{τ}_{[j,j+Nx],:} * dB^{τ}_{[j,j+Nx],:}/dθ_{j}
"""
function pbpϕ_calc_tapemap(sp::Splitting, ϕs::Vector{Float64})
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 4*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        #pm = kron([1 0; 0 -1], Ui.*Diagonal(ϕs))
        #ehk2 = exp(sp.dt*pm)*ehk2
        for xi in Base.OneTo(Nx)
            #ehkslice = ehk2[[xi, xi+Nx], :]
            ehkslice = [exp(sp.dt*Ui[xi]*ϕs[xi]) 0; 0 exp(-sp.dt*Ui[xi]*ϕs[xi])]*ehk2[[xi, xi+Nx], :]
            funcforward1 = (mu) -> bmat_Quad1ADX([exp(-sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(-2), Ui[xi], ϕs[xi], sp.dt)
            funcforward2 = (mu) -> bmat_Quad1ADX([exp(-sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(-1), Ui[xi], ϕs[xi], sp.dt)
            funcforward3 = (mu) -> bmat_Quad1ADX([exp(-sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(1), Ui[xi], ϕs[xi], sp.dt)
            funcforward4 = (mu) -> bmat_Quad1ADX([exp(-sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(2), Ui[xi], ϕs[xi], sp.dt)
            jacxi1 = derivative(funcforward1, ϕs[xi])
            jacxi2 = derivative(funcforward2, ϕs[xi])
            jacxi3 = derivative(funcforward3, ϕs[xi])
            jacxi4 = derivative(funcforward4, ϕs[xi])
            tapemap[xi, 4*si-3] = AutoJacTape([ϕs[xi]], jacxi1; val=-2, ϕ=ϕs[xi])
            tapemap[xi, 4*si-2] = AutoJacTape([ϕs[xi]], jacxi2; val=-1, ϕ=ϕs[xi])
            tapemap[xi, 4*si-1] = AutoJacTape([ϕs[xi]], jacxi3; val=1, ϕ=ϕs[xi])
            tapemap[xi, 4*si] = AutoJacTape([ϕs[xi]], jacxi4; val=2, ϕ=ϕs[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵对格点i上面的phi的导数
"""
function pbpϕ_calc_xi(sp::Splitting, hsxi::Vector{Int}, xidx::Int, ϕs::Vector{Float64}; tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpmu2 = zeros(Nt, 2, 2*Nx)
    pbipmu2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    c2t = Dict(-2 => -3, -1 => -2, 1 => -1, 2 => 0)
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = tapemap[xidx, 4*segl+c2t[hsxi[bi]]]
            jac1 = reshape(tape.output, 4, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpmu2[bi, :, :] = jac1[1:2, :]
            pbipmu2[bi, :, :] = jac1[3:4, :]
        end
    else
        @warn "计算jacobian"
        error("not implement")
    end
    return pbrpmu2, pbipmu2
end


"""
计算lnS对ϕ的导数
"""
function meas_gradϕ(ss::ScrollSVD{T}, allbmats, sp, allcfgs, ϕs;
    tapemap=missing, ignorew=false, ignorez=false) where T
    #
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    #
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    Uibar2 = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        #这里是lnw对phi的导数的B部分
        rpU, ipU = pbpϕ_calc_xi(
            sp, allcfgs[:, xi], xi, ϕs;tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            Uibar2[xi] += sum(rpU[bi, 1, :] .* real(bbar[bi, xi, :]))
            Uibar2[xi] += sum(rpU[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            Uibar2[xi] -= sum(ipU[bi, 1, :] .* imag(bbar[bi, xi, :]))
            Uibar2[xi] -= sum(ipU[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
        #
        if ignorew
            continue
        end
        #
        for bi in Base.OneTo(Nt)
            #wgt = γ*exp(sqrtma*η*ϕ)
            #注意dlnS = dlnZ - dlnw，这个负号在B矩阵时是包含在bbar中的
            #现在包含在导数的结果里
            #-∂lnwgt/∂ϕ = -(∂lnγ + ∂sqrtma*η*ϕ)/∂ϕ = -sqrtma*η
            Ui = sp.Uit[sp.vecl[bi], xi]
            sqrtma = sqrt(sp.dt*Ui/2)
            Uibar2[xi] -= sqrtma*ηdict[allcfgs[bi, xi]]
        end
        #
        if ignorez
            continue
        end
        #
        # 根据bmat_test.jl中的结果，这个相位对Z是有确实的影响的，
        # 引入phi以前Z = Z_irr
        # 引入了phi以后会变成Z_irr e^{\beta*0.5*U*phi^2}
        # d ln Z = d ln e^{-\beta (H0 - 0.5U(nu - nd + phi)^2 + Ha)}
        # d ln Z = d ln e^{-\betaH0 -\beta*(U*nud -0.5Unu-0.5Und - 0.5Uphi^2)}
        # d ln Z = d ln Z_irr e^{\beta*0.5Uphi^2 }
        # d ln Z / d \phi = d ln e^{\beta*0.5Uphi^2 } / d \phi = d \beta*0.5Uphi^2 / d\phi
        # = \beta*U*phi = ∑_{Nt} dt*U*phi
        # 这将导致phi向非0靠拢
        for bi in Base.OneTo(Nt)
            Ui = sp.Uit[sp.vecl[bi], xi]
            Uibar2[xi] += sp.dt*Ui*ϕs[xi]
            #@show sp.dt*Nt Ui ϕs[xi]
            #@show sp.dt*Nt*Ui*ϕs[xi]
        end
    end
    return Uibar2
end


"""
计算包含额外权重的导数
``````
dlnw = dlnη + dlndet(I+B) = d∑_{τ,i}ln η(τ,i) + dlndet(I+B)

dlnw / dθ_{j} = d∑_{τ,i}ln η(τ,i)/dθ_{i} + dlndet(I+B)/dθ_{j}
= d∑_{τ}ln η(τ,j)/dθ_{j} + ∑_{τ}dlndet(I+B)/dB^{τ}_{[j,j+Nx],:} * dB^{τ}_{[j,j+Nx],:}/dθ_{j}
"""
function pbpϕ2_calc_tapemap(sp::Splitting, ϕs::Vector{Float64})
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 4*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        #pm = kron([1 0; 0 -1], Ui.*Diagonal(ϕs))
        #ehk2 = exp(sp.dt*pm)*ehk2
        for xi in Base.OneTo(Nx)
            #ehkslice = ehk2[[xi, xi+Nx], :]
            ehkslice = [exp(-sp.dt*Ui[xi]*ϕs[xi]) 0; 0 exp(-sp.dt*Ui[xi]*ϕs[xi])]*ehk2[[xi, xi+Nx], :]
            funcforward1 = (mu) -> bmat_Quad2ADX([exp(sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(-2), Ui[xi], ϕs[xi], sp.dt)
            funcforward2 = (mu) -> bmat_Quad2ADX([exp(sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(-1), Ui[xi], ϕs[xi], sp.dt)
            funcforward3 = (mu) -> bmat_Quad2ADX([exp(sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(1), Ui[xi], ϕs[xi], sp.dt)
            funcforward4 = (mu) -> bmat_Quad2ADX([exp(sp.dt*Ui[xi]*mu) 0; 0 exp(sp.dt*Ui[xi]*mu)]*ehkslice, Val(2), Ui[xi], ϕs[xi], sp.dt)
            jacxi1 = derivative(funcforward1, ϕs[xi])
            jacxi2 = derivative(funcforward2, ϕs[xi])
            jacxi3 = derivative(funcforward3, ϕs[xi])
            jacxi4 = derivative(funcforward4, ϕs[xi])
            tapemap[xi, 4*si-3] = AutoJacTape([ϕs[xi]], jacxi1; val=-2, ϕ=ϕs[xi])
            tapemap[xi, 4*si-2] = AutoJacTape([ϕs[xi]], jacxi2; val=-1, ϕ=ϕs[xi])
            tapemap[xi, 4*si-1] = AutoJacTape([ϕs[xi]], jacxi3; val=1, ϕ=ϕs[xi])
            tapemap[xi, 4*si] = AutoJacTape([ϕs[xi]], jacxi4; val=2, ϕ=ϕs[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵对格点i上面的phi的导数
"""
function pbpϕ2_calc_xi(sp::Splitting, hsxi::Vector{Int}, xidx::Int, ϕs::Vector{Float64}; tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpmu2 = zeros(Nt, 2, 2*Nx)
    pbipmu2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    c2t = Dict(-2 => -3, -1 => -2, 1 => -1, 2 => 0)
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = tapemap[xidx, 4*segl+c2t[hsxi[bi]]]
            jac1 = reshape(tape.output, 4, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpmu2[bi, :, :] = jac1[1:2, :]
            pbipmu2[bi, :, :] = jac1[3:4, :]
        end
    else
        @warn "计算jacobian"
        error("not implement")
    end
    return pbrpmu2, pbipmu2
end


"""
计算lnS对ϕ的导数
"""
function meas_gradϕ2(ss::ScrollSVD{T}, allbmats, sp, allcfgs, ϕs;
    tapemap=missing, ignorew=false, ignorez=false) where T
    #
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    #
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    Uibar2 = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        #这里是lnw对phi的导数的B部分
        rpU, ipU = pbpϕ2_calc_xi(
            sp, allcfgs[:, xi], xi, ϕs;tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            Uibar2[xi] += sum(rpU[bi, 1, :] .* real(bbar[bi, xi, :]))
            Uibar2[xi] += sum(rpU[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            Uibar2[xi] -= sum(ipU[bi, 1, :] .* imag(bbar[bi, xi, :]))
            Uibar2[xi] -= sum(ipU[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
        #
        if ignorew
            continue
        end
        #
        for bi in Base.OneTo(Nt)
            #wgt = γ*exp(sqrtma*η*ϕ)
            #注意dlnS = dlnZ - dlnw，这个负号在B矩阵时是包含在bbar中的
            #现在包含在导数的结果里
            #-∂lnwgt/∂ϕ = -(∂lnγ + ∂sqrtma*η*ϕ)/∂ϕ = -sqrtma*η
            #-∂lnabs(wgt)/∂ϕ = -∂ln abs(e^{(rsq+im*isq)*η*ϕ}) /∂ϕ = -∂lne^{rsq*η*ϕ}/∂ϕ
            # =-rsqma*η
            Ui = sp.Uit[sp.vecl[bi], xi]
            sqrtma = sqrt(complex(-sp.dt*Ui/2, 0.0))
            Uibar2[xi] -= real(sqrtma)*ηdict[allcfgs[bi, xi]]
        end
        #
        if ignorez
            continue
        end
        # 根据bmat_test.jl中的结果，这个相位对Z是有确实的影响的，
        # 引入phi以前
        # 为Z = Z_irr e^{-\beta*0.5*U}，也就是e^{\beta*0.5*U}*Z = Z_irr，
        # 也就是bmat_test.jl中乘exp(0.4*1.5/2)的原因
        # 引入了phi以后会变成Z_irr e^{-\beta*0.5*U} e^{-\beta*0.5*U*phi^2}
        # d ln Z = d ln e^{-\beta (H0 + 0.5U(nu + nd - 1 + phi)^2 + Ha)}
        # d ln Z = d ln e^{-\betaH0 - \beta*(U*nud -0.5Unu-0.5Und + 0.5Uphi^2 - Uphi + 0.5U)}
        # d ln Z = d ln Z_irr e^{-\beta*0.5*U} e^{-\beta*0.5Uphi^2 + beta*Uphi }
        # d ln Z / d \phi = d ln e^{-\beta*0.5Uphi^2 + beta*Uphi } / d \phi = d (-\beta*0.5Uphi^2 + beta*Uphi) / d\phi
        # = -\beta*U*phi + beta*U = ∑_{Nt} -dt*U*(phi-1)
        # 这将导致phi向1靠拢
        for bi in Base.OneTo(Nt)
            Ui = sp.Uit[sp.vecl[bi], xi]
            Uibar2[xi] += -sp.dt*Ui*(ϕs[xi]-1)
        end
    end
    return Uibar2
end

